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Abstract. In an introductory course on dynamical systems or Hamiltonian mechanics, vector 
field diagrams are a central tool to show a system’s qualitative behaviour in a certain domain. 
Because of their low sampling rates and the involved issues of vector normalization, these 
plots give only a coarse insight and are unable to convey the vector field behaviour at 
locations with high variation, in particular in the neighborhood of critical points. Similarly, 
automatic generation of phase portraits based on traditional sampling cannot precisely capture 
separatrices or limit cylces. In this paper, we present ASysViewer, an application for the 
interactive visual exploration of two-dimensional autonomous dynamical systems using line 
integral convolution techniques for visualization, and grid-based techniques to extract critical 
points and separatrices. ASysViewer is addressed to undergraduate students during their first 
course in dynamical systems or Hamiltonian mechanics. 

(Some figures may appear in colour only in the online journal) 


PACS numbers: 01.50.H-,01.50.hv 


1. Introduction 

The qualitative behaviour of two-dimensional autonomous dynamical systems is traditionally 
explored visually by means of either a vector held diagram or as a phase portrait that shows 
critical points, several typical orbits, and separatrices. Both techniques, however, have the 
disadvantage that they only coarsely sample the region of interest. Furthermore, traditional 
automatic generation of these plots, in general, cannot precisely capture the separatrices or 
reproduce limit cylces. A continuous visualization of the phase space can be achieved using 
line integral convolution (LIC), which is a well-known technique for vector held visualization 
introduced by Cabral and Leedom III in 1993. Unfortunately, up to now, this technique has 
not found its way into lectures about dynamical systems because of the lack of software that 
can produce high-quality phase images interactively. 

The aim of this article is to present how visualization and numerical techniques can help 
in the visual and interactive exploration of two-dimensional autonomous dynamical systems 
which might be helpful in undergraduate courses on classical Hamiltonian mechanics as 
well as on dynamical systems in general. For that, we mainly use line integral convolution 
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techniques, integral curves, slicing, and topology extraction, which we bring to the graphics 
board. 

A detailed introduction to general dynamical systems can be found in the standard 
literature, e.g.. Lynch J2| or Wiggins 0. A state of the art report in flow visualization 
using dense and texture-based techniques, in particular various LIC techniques, was given 
by Laramee et al. 0. 

The structure of this paper is as follows. In section[2] we briefly review 2D autonomous 
systems and the classification of critical points. In section [3] we discuss the standard 
vector field diagram and how it can be replaced by line integral convolution visualization. 
Furthermore, we detail our topology extraction method. In section [4] we present our 
application and give several examples in section[5] 

Our software ASysViewer is based on the cross-platform framework Qt |5| and the Open 
Graphics Library (OpenGL) (6). The source code and Windows binaries are available from 
http://go.visus.uni-Stuttgart.de/asysviewer 


2. Autonomous dynamical systems 


A two-dimensional autonomous dynamical system is defined by two coupled first-order 
differential equations 

x = P(x,y), y = g(x,y), (1) 


where a dot represents differentiation with respect to time, f, and P and Q are two (arbitrary) 
C 1 -differentiable functions of x and y in an open subset G C R 2 . Both, x and y themselves 
are functions of t. Because of the theorem of existence and uniqueness, solutions of Eq. 0 
do not cross. Hence, we can interpret the system of equations (|T]i geometrically as a two- 
dimensional vector field with a vector /(x,y ) = (P(x,y),Q(x,y)) T attached at each point 
(x,y) € G. The trace of a solution (7(t,xo,yo) = (^x(t,xo,yo),<y y (t,xo,yo)) T of the coupled 
system 0 with initial values (xq ,yo) and parameter t is called integral curve (trajectory, orbit, 
flow line, characteristic line, or streamline), where <7 is always tangential to the vector / at 
the current point (x,y). Contour lines for which Ay/Ax = y/x = const are called isoclines. 

To obtain a qualitative view of the vector field, a phase portrait that consists of critical 
points and several characteristic lines and separatrices is used. Stationary points (also 
called fixed or equilibrium points) are defined by x = y = 0. Critical points are isolated 
stationary points, i.e., stationary points where the determinant of the Jacobian does not vanish, 
det(/) / 0. The behaviour of the system 0 in a close neighborhood of a critical point (x c ,y c ) 
is characterized by the eigenvalues X \2 of the Jacobian 
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and X = x —x e , Y y -y c . Then, the eigenvalues A 1.2 can be determined by means of the 
characteristic equation det(7 AI 2 ) = 0, 


^ 1,2 = 


tr(/)± V /(tr(/)) 2 -4det(7) 


(4) 


with trace, tr(7) = J\ \ +J22, determinant det(7) = J11J22 / 12/21 ^ and I 2 being the identity 

matrix in two dimensions. If the real part of the eigenvalues is nonzero, the critical point 
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is called hyperbolic (structurally stable, i.e., robust against perturbations), and, according 
to Hartman’s theorem, the phase portrait of the original system in the neighborhood of this 
critical point resembles that of the linearized system. Figure [T] shows the classification of 
critical points depending on the values of tr(7) and det(7). If det(7) < 0, the critical point can 
only be a saddle point. If det(7) > 0, it depends on tr(7) whether the critical point is stable or 
unstable, i.e., in which way orbits approach or move away from the critical point, respectively. 



Figure 1 . Classification of critical points depending on the trace and the determinant of the 
Jacobian J (see also Lynch(2))- The parabola is given by det(7) = (tr(7)) 2 /4, see Eq. jJJ. 


Figure [2] shows examples of the six main types of critical points which are located at the 
centers of the vector field diagrams. The arrows are normalized to unit length to better see the 
characteristics of the field. The red curves represent some typical orbits. 

To further characterize a vector field, one has to determine the characteristic lines that 
mark the boundary between different regions. These lines are also called separatrices, in 
general. A characteristic line is called homoclinic orbit if the limits 

lim a(/,xo) = lim d(t,xo)=x c (5) 

have the same critical point x c . If the limits have different critical points, the characteristic 
line is called heteroclinic orbit', otherwise, it is a generic orbit. 

3. Visualization techniques 

3.1. Vector field diagram 

The most straightforward visualization of the two-dimensional system 0 is a vector field 
diagram with arrows pointing in the direction / = ( P,Q) T as already shown in figure [ 2 ] 
Figures |3(a)|(b)| show a vector field diagram where the arrows are either normalized or not. 
Without normalization (Fig. |3(a)| , the arrow lengths are too diverse so that the overall structure 
of the vector field g ets los t if the diagram cannot be scaled interactively. 

While in figure [3(b)] the critical points at ?i = (— 1,0 ) 7 and C 2 = (0,0) r are obviously a 
center and a saddle, respectively, the location and the types of the critical points in figure [3(c)] 
C 3 = (0,0) r , C 4 = (0.734, —1.363) r , C 5 = (1.931,— 0.518) r , are not as obvious. Even more 
complicate is the retrieval of the separatrices in figure [3(c)| 

3.2. Line integral convolution 

Line integral convolution is a well-known texture synthesis technique in flow visualization. It 
has the advantage that its spatial resolution is much higher than that of a vector field diagram. 
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(a) unstable node 


(b) stable node 


(c) saddle 




(d) unstable focus 


(e) center 


(f) stable focus 


Figure 2. Vector plots for critical points at ( x c — 0,y c = 0). The arrows are normalized to unit 
length. The red lines show some example orbits. 



(a) (b) (c) 


Figure 3. Vector field diagram for P(x,y) = y, Q(x,y) = x + x 2 with critical points c\ = 
(—1,0)' and ?2 = (0,0)'; (a) real length, (b) normalized, (c) Normalized vector field 
diagram for P(x,y) — Ax — 2x 2 — y 2 , Q(x,y) — x(l +xy) with critical points C 3 = (0,0) r , 
?4 = (0.734, —1.363) r , andcs = (1.931, —0.518) . 


The idea is to locally smear an image (usually an image consisting of white noise) along the 
direction of a vector field. Thereby, pixel values along streamlines are strongly correlated 
while orthogonal to the vector field, there is virtually no correlation. 

The basic algorithm works as follows. First, we have to generate a window-filling noise 
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texture T of size N x x N y pixels with random numbers in the range [0,1] for each pixel that 
covers the domain £2, see figure [4] 



(a) (b) 


Figure 4. (a) 2D domain of interest Q = {(x min ,y min ), (x max ,;ym a x)}- (b) Noise texture/image 
T with resolution N x x N y pixels, slightly smoothed by a Gaussian kernel to reduce high 
frequencies. 

Next, we perform a discrete line integral convolution of the noise image for each pixel 
p G T using the trajectory d(t,xo). The resulting intensity I (p) at the pixel p follows from 

N 

I( p) = *(0)r(p) + £ [k(n)T(u n ) +k(-n)T(y n )] , (6) 

n— 1 

where A: is a filter kernel and 2N + 1 is the filter length. In detail, we start with the intensity 
T (p) £ [0,1] of the noise texture at pixel p weighted by k( 0) and then sum up all the intensity 
values of T along the positive and negative parts of the trajectory weighted by the kernel value 
k at the corresponding sampling points u„ and v„, respectively.” For that, we first need the 
relation between pixel coordinates q = (q x ,q y ) £ T and domain coordinates x = (x.y) 1 £ £2 
which is defined by 

x i->- q with q x = N x -, q y = N y -. (7) 

-^max -^min ymax ymin 

For the sake of convenience, we employ the Euler method to integrate the trajectories 
= vo) i n positive and negative direction. The initial pixel p defines the initial positions 
«o = vo which follow from the inverse of equation 0- Then, Euler iteration delivers the 
sampling points u n ,v n £ £2 with 

m„+i =u n +hf n , v„+i =v„-hf„, n = 0,...,N—l. (8) 

According to equation 0. we obtain the sampling points u„ , v„ £ T from u n and v„. 

The filter kernel for standard LIC can be any symmetric function. Here, we use a Hann 

filter 

1 / 7 Tfi \ f- 

2NV +C0S ~n) ’ n = E Kn) = \. (9) 

The line integral convolution is computationally quite expensive. But with the high 
parallelism of GPUs, we can accomplish the convolution at interactive rates. At the end, 
we obtain the LIC image composed of all the intensity values I (p) calculated above, which 
can afterwards be scaled arbitrarily to enhance the field structures, see, e.g., figure[6] 

Unfortunately, the standard LIC loses the directional information of the vector field, i.e., 
forward and reverse direction are indiscernible. To resolve this disadvantage, we can use an 
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asymmetric filter kernel that comprises orientation (oriented LIC). see, e.g., Wegenkittl et 
al. Q. The most straightforward kernel for that is a sawtooth. Furthermore, for oriented LIC, 
we have to replace the noise texture of figure[4]by a texture with a finite number of randomly 
distributed spots and black background. Figure [5] shows the vector fields of figure [3] using 
oriented LIC. Additionally, we use a repeating colour map for the vector magnitudes ||/||. 



Figure 5. Oriented line integral convolution images of the vector fields of figure [3] generated 
by ASvsViewer using the scripts paperExpl. js and paperExp2. js. Here, colour encodes 
vector magnitude ||/||. As these values are very diverse, we repeat the colour map with n 6 N. 


Another possibility to maintain the direction information of the vector field is to animate 
the line integral convolution (animated LIC). For that, we change the frequency of the Hann 
filter and modulate the phase by animation time T. Additionally, we multiply the modified 
Hann filter with a Hann window. Thus, k(n) from equation (|9|) is replaced by 


h(n) 


1 

4 


1 T cos 



(l+0.2 cos^). 


( 10 ) 


3.3. Topology extraction 

So far, we have demonstrated the benefits of dense vector field visualizations by the LIC 
family of algorithms as compared to the discrete sparse visualizations by vector plots: due 
to their dense nature, they provide detailed information about the direction of /, and to some 
extent also about its magnitude and orientation. However, it is still difficult to spot the critical 
points and, in particular, the separatrices emanating therefrom (figure |5j. In other words, it is 
still difficult to identify regions of qualitatively different behaviour. 

Topology extraction usually starts with the extraction of critical points. A major difficulty 
with this step is, however, the lack of a theory on the existence and location of critical points 
in arbitrary (higher-order) vector fields. Nevertheless, in the special class of low-order vector 
fields arising from 2D (bilinear) tensor-product linear interpolation, their number is restricted 
to two. Hence, in vector fields where Cl is discretized into cells and / is discretized at the nodes 
(the corners of the cells), such a field is present in each cell due to bilinear interpolation. 

The algorithm to detect the critical points works as follows: 

(i) Detect cells with different signs at their nodes both in terms of P and Q. 







Visual exploration of 2D autonomous dynamical systems 


7 


(ii) Such cells are candidates for stationary points. However, a sationary point is present only 
if the zero-level isocontours of P and Q intersect. 

(iii) Analytical derivation of these intersections is possible, but subject to numerical 
instability. Thus, we follow a subdivision approach: 

(iv) Each cell with different signs in P and Q is subdivided, and the resulting four children 
cells are again tested for different signs in P and Q. This procedure is repeated until no 
cell exhibits different signs in P and Q or until a maximum subdivision level (a prescribed 
accuracy) is achieved. The centers of the remaining cells with different signs represent 
stationary points. 

(v) Those stationary points that have det(J) f 0 are the desired critical points. 

Those critical points where det(7) < 0 are of type saddle and hence give rise to 
separatrices. To escape the zero-velocity region of the critical point, for each of the four 
separatrices an offset step along the major (or minor) eigenvector of the Jacobian is applied 
and streamline integration in forward (or reverse) direction is carried out, resulting in the 
desired separatrices. 


4. The viewer 


4.1. Graphical user interface 

A screenshot of ASysViewers graphical user interface is given in figure [6] The main window 
shows standard LIC with critical points and separatrices of the “mathematical pendulum with 
friction” example, see Eq. (111. In the “System” window, the functions P(x,y) and Q(x.y) 


can be given as mathematical expressions, and the domain 12 can be controlled either using 
mouse interaction (zooming: middle button; panning: left button) or by typing the values 
in the corresponding boxes. With the right mouse button a trajectory can be shown which 
starts at the current mouse position. The parameters for this trajectory can be set within the 
“Trajectory” window. The “FuncPlot” window shows the values of P and Q along the 
white intersection line drawn in the main window. This line can be set with the right mouse 
button and “intersection” selected in the “Show” window. After defining the dynamical 
system, the critical points as well as the separatrices are determined by pressing the “calc” 
button in the “Topology” window. 

There are also two windows to control the parameters a and eps which can be used as free 
parameters to define the functions P and Q. 


4.2. Technical details 

ASysViewer is based on the cross-platform application and user interface framework Qt 0, 
the graphics library OpenGL, and the OpenGL Shading Language (GLSL) |6|. In the main 
window, we draw a single quad that covers the whole 2D domain of interest, see also figure [4] 
For each fragment (pixel) of this quad, the line integral convolution is calculated within a 
fragment program/shader that works as a SIMD (single instruction multiple data) architecture. 
The fragment program is written in a C-like programming language (GLSL) and is compiled 
and transferred to the GPU at runtime. When the functions P and Q are modified, the function 
strings are inserted in the fragment code and the fragment program is recompiled. Due 
to limitations of mathematical functions/operators in GLSL, expressions like “x 4 ” must be 
formulated as “x *- x * x”. 
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Figure 6. Screenshot of ASysViewer 's graphical user interface showing standard LIC with 
separatrices and critical points. In the “FuncPlot” window, the blue and red lines correspond 
to the functions P and Q evaluated along the interactively defined white line in the main 
window. 


The integration of the trajectories to visualize the separatrices or individual orbits 
selected via mouse control is done using a step-controlled Runge-Kutta-Fehlberg integrator. 
For that, the functions P and Q are parsed using the muparser ( 8 ) library by Ingo Berg. For 
better reproducibility, the parameter settings of ASysViewer are also scriptable by means of 
QtScript which is based on the ECMAScript J9] standard. The necessary scripting functions 
can be checked up from the source code. 

5. Examples 

5.1. Mathematical pendulum with friction 

As a first example, we consider the damped oscillation of the mathematical pendulum of unit 
length which is described by the equation of motion (p + ft) 2 sintp + atp = 0 with deflection 
angle (p. angular frequency to, and a positive frictional coefficient a > 0 , see figure [ 6 ] 
Substituting (p = x, <p =y converts the second order ordinary differential equation into the 
standard form of equation 

x = P(x,y) =y, y = Q{x,y) = —orsinx — ay. (11) 

The positions of the corresponding critical points immediately follow from x = 0 = y. Thus, 
y c = 0 and x c = sin(H 7 f) with n £ TL. The Jacobian yields tr(J) = a and det(7) = ft ) 2 cos(« 7 r). 
If n is odd, the determinant is negative and the critical point is a saddle. If n is even, on the 
other hand, the critical point is either a center (undamped oscillator, a = 0 ), an unstable focus 
(a 2 < 4 ft) 2 cos(« 7 r)), or an unstable node. In the damped case, each unstable focus/node 
is connected to the two neighbouring saddle points via heteroclinic orbits, whereas in the 
undamped case, each saddle point is connected to its two neighbouring saddle points. 
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A limit cylce is an isolated periodic solution which is defined by the set of all points y with 
o(t , x) —» y for t —> <=°. An example can be found in the phase space for the differential equation 
related to the oscillation of a violin string as derived by Rayleigh Qo), 


v-p 



v + v = 0, 


P >0, 


( 12 ) 


which can be brought into the standard form (Eq. (|T]i) via v = x and v = y. 
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Figure 7. (a) Phase space image of the Rayleigh system, Eq. GU. visualized using oriented 
LIC. The trajectory (yellow line) starting at x = 0.5,v — 0.2 approaches the limit cycle 
asymptotically. Here, Q. = [—2.5,2.5] X [—2.5,2.5] and p = 1. (b) Vector field diagram. 


As shown in the phase space image, figure 7(a) the oriented LIC representation of the 
Rayleigh system indicates that there must exist a boundary curve separating the outflowing 
and inflowing vector field. Every integral curve a(t,x /) approaches this boundary curve (limit 
cycle) irrespective of the initial point 3c 


5.3. Poincare index 


The Poincare index ip of a closed curve T is determined by integrating the change in the angle 
of the vectors at each point of F along T. For the numerical calculation, we set F to a circle 
of radius r and replace the integration by the sum 


ip — 


1 

2n 


N 

^ arcsin 

n =0 


II fn x fn+l |1 

WfnUfn+lW’ 


(13) 


where f n = f (r cos %,r sin (p n ) T and % = 2nn/N. For a sink, source, or center, we have 
ip = +1. A saddle point has ip = — 1. For details to the Poincare index, we refer the reader 
to Wiggins j3j. In ASysViewer. the Poincare index of an arbitrary point x. can be determined 
interactively by constructing a circle around x via mouse handling (right button). For that, 
select “poincare idx” in the mouse combo box of the “Show” window. 
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5.4. Pitchfork bifurcation 
Consider the system 

x = ax— x 3 , y = -y, (14) 

with free parameter a £ R. If a is continuously varied from negative to positive values, this 
system undergoes a pitchfork bifurcation. As long as a is negative, there is only a single 
critical point at x neg = (0.0) 7 . For a = 0, there is one non-hyperbolic critical point at the 
origin. If a > 0, we obtain three critical points x\ = (0,0) r (saddle), xi = (—-y/ a,0) T (stable 
node), and j ?3 = (i/a, 0 ) r (stable node). 

Figure[ 8 ]shows the “pos,neg P. Q” image for the pitchfork bifurcation system with a = 1, 
where colour indicates the signs of the functions P and Q. Hence, a point where all four 
colours hit is a stationary point, P = Q = 0. 




Figure 8. The system (Eq. 114) ) shows a pitchfork bifurcation at a = 0. If a > 0, the critical 
point at the origin is a saddle whereas the other two are stable nodes. The colour coding in 
(a) is as follows: yellow (P > 0,Q > 0), red (P > 0,Q < 0), brown (P < 0, Q < 0), green 
(P < 0, Q > 0); (b) pitchfork bifurcation. 


6. Summary 

In this work we presented a framework for interactive visual exploration of 2D autonomous 
dynamical systems. For that, we use different types of line integral convolution techniques for 
a dense representation of the vector field. We also described a robust algorithm to numerically 
detect critical points. As the source code of ASysViewer is freely available, the interested 
student is encouraged to conduct his or her own experiments and to look also at the coding 
details. 
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